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We consider the simplest generalizations ol the valence bond physics of SU (2) singlets to SU (N) 
singlets that comprise objects with N sites — these are SU(N) singlet plaquettes with N = 3 and 
N = 4 in three spatial dimensions. Specifically, we search for a quantum mechanical liquid of such 
objects — a resonating singlet valence plaquette phase that generalizes the celebrated resonating va- 
lence bond phase for SU(2) spins. We extend the Rokhsar-Kivelson construction of the quantum 
dimer model to the simplest SU (4) model for valence plaquette dynamics on a cubic lattice. The 
phase diagram of the resulting quantum plaquette model is analyzed both analytically and numer- 
ically. We find that the ground state is solid everywhere, including at the Rokhsar-Kivelson point 
where the ground state is an equal amplitude sum. By contrast, the equal amplitude sum of SU(3) 
singlet triangular plaquettes on the face centered cubic lattice is liquid and thus a candidate for 
describing a resonating single valence plaquette phase, given a suitably defined local Hamiltonian. 



I. INTRODUCTION 

The physics of valence bondsi^ has provided a fruitful 
framework for understanding non-Neel phases of SU (2) 
invariant quantum magnets in low dimensions. 3 "^ In par- 
ticular, its reformulation as the study of quantum dimer 
models by Rokhsar and Kivelson 5 - has yielded a sharp 
but intuitive description of a diverse set of unconven- 
tional phenomena, most notably the existence of topo- 
logically ordered, fractionalised RVB liquids 6 -! 7 - originally 
proposed by Anderson and the phenomenon of Cantor 
deconfmemenl£. Extensions of such dimer models to 
higher dimensions have been considered 9 - ! 10 ' 11 ! 12 ' 13 ! 14 and 
in the process the existence of two distinct RVB liquids in 
d = 3 has been demonstrated. The dimer model formu- 
lation has the additional virtue of having the manifest 
interpretation of a strongly coupled, frustrated, gauge 
theory with RVB phases arising as deconfmed phases in 
such gauge theories! 15 ! 16 It is also useful to note the re- 
cent construction of finite range spin models that realize 
dimer models arbitrarily accurately in their low energy 
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We should note that the interest in RVB liq- 
uids also comes from a much larger study of topological 
phases in condensed matter systems^ 

In this paper we pursue a generalization of the valence 
bond idea to SU(N) groups with N > 2 with the aim 
of searching for topological phases in models with these 
higher symmetries. The valence bond is replaced by a 
singlet formed by N spins in the fundamental represen- 
tation of SU(N). The basic geometrical degrees of free- 
dom thus become triangles for N — 3, and objects con- 
taining four sites, such as square plaquettes for N = 4, 
and so on. We restrict our study to singlet plaquettes 
with N = 3, 4. With the assumption that for suitable 



microscopic Hamiltonians the low energy sector is well 
described by coverings of the lattice by such objects, we 
are led to study quantum plaquette models. 

These plaquette models are similar in spirit to two- 
form gauge theories, where degrees of freedom are also 
defined on plaquettes. 20 However, there is one important 
difference. In the latter, there is a local constraint on pla- 
quettes that share a bond. In our problem, the site origin 
of the plaquette variables is reflected in their satisfying 
a site constraint — specifically, that one and only one of 
the plaquettes that share a common site be occupied by 
a singlet. 

Our strategy in this paper, mirroring that used in the 
study of valence bond physics, is to begin with a Hilbcrt 
space spanned by orthogonal states \a) labelled by ad- 
missible plaquette coverings a of the lattice under study. 
In this space we study the phase diagram of the simplest 
local Hamiltonian which includes the quantum dynamics 
of resonance between pairs of plaquettes and a potential 
energy of interaction between them which we shall term 
the quantum plaquette model (QPM). As is well known, 
when the signs of the matrix elements are of the "Perron- 
Frobenius" form (i.e. no positive off-diagonal matrix el- 
ements), such models exhibit a Rokhsar-Kivelson point 
(which occurs at v — 2t for the case of the cubic lattice, 
Eq. 0} where the equal amplitude superposition, 



(i) 



saturates a lower bound on the ground state energy of 
the problem. The diagonal correlations in this state de- 
fine a classical statistical mechanics of plaquettes which 
can be studied to establish the true nature of this poten- 
tially liquid state. Knowledge of this special point, along 
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with the analysis of other limits typically allows the phase 
diagram to be established with some confidence. The re- 
maining challenge is to show that one can actually define 
a sign convention for wavefunctions of the spins that lead 
to the desired signs of the matrix elements and finally to 
construct local spin Hamiltonians that encode all of this 
physics. 

For the N = 4 problem on the cubic lattice we report 
considerable progress on this strategy. We show that 
matrix elements of the nearest neighbor spin interaction 
have the desired form. We analyze the statistics of the 
RK point wavefunction computationally, using a plaque- 
tte version of the pocket Monte Carlo algorithm devel- 
oped for dimers in Ref. [U We find that the classical 
system exhibits long range order, in which the prefer- 
able plaquette locations form a cubic lattice, analogous 
to the plaquette valence bond crystal for dimers. This 
crystalline phase is accounted for analytically from sev- 
eral perspectives. Besides this plaquette phase, the phase 
diagram of the QPM also contains the staggered and 
columnar phases and thus the model does not allow a 
topological phase. By contrast, we find that the equal 
amplitude sum for the N — 3 problem on the face cen- 
tered cubic (fee) lattice is liquid and thus a candidate 
for an RSVP phase. The caveat here arises from the 
non-ergodic nature of the dynamics in the simplest QPM 
on this lattice. While we suspect that a longer ranged 
(but still local) dynamics will fix this, it remains an open 
problem at this time. 

We note that we are, by no means, the first to con- 
sider such higher symmetry problems and the emergence 
of plaquette variables — the new element in our work is 
the study of QPMs and the idea of the plaquette liq- 
uid. In d — 1, SU(4) symmetric models are readily 
constructed by fine-tuning an SU(2)xSU(2) spin-orbital 
models , 22 ! 23 ' 2 or a two-leg ladder.— In two dimensions, 
singlet plaquettes have already been considered on a tri- 
angular lattice i 22 ' 26 Separately, SU(4) models have been 
studied en route to connecting large- N SU(N) models 
with the SU(2) case i 27 ' 28 The prospect of studying spin- 
3/2 cold atom systems has inspired an elegant body 
of work that includes relevant results on SU(4) pla- 
quette solids and even some generalizations to SU(N) 
ladder s 29 ' 30 . It is also useful to remark here that models 
with interspersed conjugate representations on bipartite 
lattices^ exhibit singlet dimers ( "mesons" ) while we con- 
sider models with uniform choices of representation which 
are thus forced to exhibit larger objects ("hadrons"). 

This paper is organized as follows. In Sect. Ill Al we 
repeat Rokhsar and Kivelson's original "derivation" of 
the N = 4 QPM — essentially, it shows how the appro- 
priate signs can be generated for the matrix elements of 
the QPM. In Sect. Ill Bl we briefly introduce the pocket 
algorithm and present the results of numerical simula- 
tion supporting our view of the ground state at the RK 
point of the model. In Sect. Ill CI we provide an analyt- 
ical explanation of the plaquette phase based on a low- 
and high-dimensional generalisations of the model. In 
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FIG. 1: Minimal phase diagram of the singlet valence plaque- 
tte model on a cubic lattice. The transitions are first order. 

Sect. IIIII we present an analysis of the N — 3 QPM on 
the face-centred cubic lattice. In the conclusion, we com- 
ment on the model's utility and open questions. 

II. THE N=4 QPM ON THE CUBIC LATTICE 
A. Microscopies and Hamiltonian 

The quantum plaquette model is an effective model 
for the SU(4) Heisenberg antiferromagnet in a valence 
plaquette dominated phase, in much the same way the 
quantum dimer model is an effective model for the SU(2) 
Heisenberg antiferromagnet in a valence bond dominated 
phase. As we have noted above, the choice of plaquettes, 
rather than dimers, as new degrees of freedom is dictated 
by the fact that in the SU(4) model the lowest energy 
configuration on a four site lattice is a singlet which mixes 
equally all four sites. 

In this section we will establish some useful properties 
of the microscopic wavefunctions that correspond to pla- 
quette coverings of the cubic lattice. Should it prove pos- 
sible to construct a Hamiltonian in the full SU(iV) invari- 
ant Hilbert space which has the desired coverings, and no 
other states, as ground states — a non-trivial problem — 
the results in this section can be used to imitate Refs. 
[l7llT8l and rigorously realize the QPM that we introduce 
here. We will not have much to say about this possibility 
beyond a few closing remarks in the last section. 

A classical plaquette covering of the lattice corresponds 
to a quantum state, the wave function \a) of which is 
the product of individual plaquette singlet wave func- 
tions. These states are not orthogonal and computing 
the Hamiltonian matrix elements requires knowledge of 
the overlap matrix elements. In the case of the dimer 
model, a simple expression for the overlap matrix ele- 
ments was given in terms of the length of loops in the 
transition graph (which is formed by superposing the two 
dimer coverings under consideration)^ 

For the plaquette model, the situation is more involved 
as the geometry of the transition graph elements is much 
more complex than the loops formed by the dimers. We 
have depicted the simplest members in Fig. [51 and it is 
easy to check that their contribution to the overlap is 
given by a factor 

S = 6-^ +1 (2) 

Here I is the total number of plaquettes in the transition 
graph element. With this definition of I, the overlap in 
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FIG. 2: (color online) Transition graphs with various number 
of strongly connected parts. 

the SU(2) dimer case is given by 2~i +1 , which shows 
that the value of the convergence parameter, x = 1/6, 
is much smaller in the SU(4) case than it is for SU(2), 
where it equals 1/2. 

In general, configurations contributing to the overlap 
can be mapped bijectively onto an edge four-coloring 
model on the graph in which the centres of the plaquettes 
define the sites and their overlapping corners the edges. 
For the case of a transition graph element defined by two 
overlapping checkerboards, one thus finds that the over- 
lap matrix element is even lower, as that entropy^ im- 
plies .116'/ 2 . Interestingly, the Pauling estimate for this 
quantity is the same as for the ice model on the same 
graph, namely 0.094'/ 2 , although it is obvious that the 
entropy of the ice model is lower. 

We are not aware of a good rigorous bound for the 
entropy of such loop models; it is straightforward to show 

L_ 

(AppendixEJ) that S < c (si /2j 2 = 0.711^, (where c is 

an unimportant constant) is such a bound, although not 
at all a strict one. Indeed, all the cases we have explicitly 
checked have Eq. [5] as an upper bound: S < 6~2 +1 . 

It is convenient to work in the orthogonal basis \a) 
which can be constructed from the original basis |ev) as 
\a) = (5~2) QQ , |q/), where S aa i = (a\a'}. Such an or- 
thogonalization can be conveniently carried out analyti- 
cally assuming the off-diagonal overlap matrix elements 
to be small. 

For example, consider the projection of the nearest 
neighbor SU(4:) Hamiltonian 

ff = EE s :( i ) s ™W' ( 3 ) 

(ij) run 




FIG. 3: (color online) Transition graph generated by two con- 
figurations (denoted light and dark, respectively) of three pla- 
quettes. A strong connection involves two plaquettes sharing 
an edge, whereas a weak connection is formed by two plaque- 
ttes coinciding at a single site only. 

where the are SU(A) generators, to the valence pla- 
quette subspace. To this end one expands S aa i or 
Ha&> = (ev|-ff|<S') in powers of x. Retaining the lead- 
ing off-diagonal (kinetic) and diagonal (potential) terms, 
one obtains the QPM Hamiltonian (also referred to as 
the RK model): 

H = -tY,W0\+vY^\a){a\ (4) 

where t = 5x and v = 16x 2 (see Appendix [Cl for details). 

The precise numerical values for t and v should not 
be taken seriously - due to the presence of higher-order 
corrections and the absence of other simultaneous per- 
turbations that one can consider, and the choice of pla- 
quette coverings itself remains to be energetically forced. 
Nonetheless the negative sign of the kinetic term sug- 
gests that it is reasonable to study a QPM with the above 
"Perron- Frobenius" form which is what we will do next. 
For such models one has the highly useful feature that 
the ground state is a nodeless wavefunction. 

We note that the states a and (3 are connected by a 
single flip of a pair of plaquettes. It is easy to show that 
one ground state at the RK point (v = 2t) is the equal 
amplitude superposition of all possible states \a). Study- 
ing its properties reduces to sampling (with equal prob- 
ability) different plaquette coverings. We use a pocket 
algorithm 2 ^ to implement this task. 



B. Pocket algorithm and MC simulation results 

The main advantage of the pocket algorithm is that 
it efficiently finds possible cluster moves. The cluster 
moves, starting from a covering A, are constructed as 



follows. 1) Consider a transition graph, C, between A 
and another covering, B. 2) Pick one connected part 
of the transition graph, V € C; V is the content of the 
pocket. 3) Replace the plaquettes of A belonging to V 
by those of B. The new covering A' obtained in this way 
from A is an allowed covering, because it automatically 
satisfies the close packing constraint. 

Finding an uncorrelated new covering B, is of course in 
general not an easy problem. In the pocket algorithm, it 
is obtained by acting on A with a symmetry operation 34 
of the lattice (reflection, for example). Replacing A with 
a reflected version of itself would of course not yield an 
independent configuration, but partial replacements (us- 
ing V only) eventually do. 

This algorithm has been applied to the dimer 
problem^ In that case the pocket contains a one dimen- 
sional object which was shown not to invade space, as 
loops in dimer transition graphs never branch. The sit- 
uation is different with plaquettes. It turns out that the 
transition graph actually percolates and the system is not 
very close to the percolation threshold. In other words 
the typical size of non-percolating pocket is small. This, 
however, does not pose a problem, as the non-percolating 
pockets contain about 40% of all plaquettes, that is only 
about 60% of the plaquettes are flipped when the perco- 
lating pocket is encountered. 

The limiting factor in our case is the quick growth of 
the equilibration time r with increasing size of the lattice. 
The growth appears to be faster than L 3 with L being 
the size of the lattice. To establish r we used a binning 
procedure (see Appendix [D] for details) . Due to these 
ergodicity restrictions we were limited to a linear size 
L — 64, i.e. a total of 64 3 /4 = 65536 plaquettes, which 
is sufficient to establish the ordering properties of the 
model. 

In the simulations we measure various (equal time) 
correlation functions, as a function of plaquette sep- 
aration. The coordinates r = (x,y,z) are assigned 
to a plaquette covering sites (i,j,k,l) according to the 



rule: x = mm(x l ,x j ,Xk,xi), y 



nin(yi,yj,yk,yi) and 
z = ni\a.{zi, Zj, Z}~, zi), where (xi,yi,Zi) are the coordi- 
nates of the site i etc. In the same way we can assign 
coordinates to a flippable pair of plaquettes (which we 
also call a resonating cube). We use letters X, Y and 
Z to denote the orientation of the plaquette or the pair 
forming a 'resonating cube'. A letter corresponds to the 
axis orthogonal to a given plaquette (or plaquettes, in 
case of a resonating cube). Most generally we are inter- 
ested in studying the correlation function G na f){r), where 
n = 1 for the plaquette correlator, n = 2 for the resonat- 
ing cube correlator, and a and (3 denote orientations. In 
fact we will only need to consider the case of a = (3. Us- 
ing cubic symmetry we drop the orientation index and 
denote the correlation function as G n (r±,rt\), where r± 
and ry refer to the components of r orthogonal and par- 
allel to the plaquette. 

In Fig. [?1 we plot the plaquette correlation function 
Gi(0,r). It displays oscillations persisting to large sepa- 
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FIG. 4: Correlation function Gi(0,r) for L = 64. 
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FIG. 5: Correlation function Gi(r,0) for L = 64. 



ration r, so that one may conclude that the system pos- 
sesses long range order. The simplest possibility might 
be some kind of orientational order (rotational symme- 
try breaking). Quite a different conclusion may be drawn 
from the plot of Gi(r, 0), shown in the Fig. [5j At face 
value, the plot suggests an absence of the long range or- 
der, because the oscillations of the correlation functions 
decrease fast with distance. To detect a possible pres- 
ence of orientational order we have computed the quan- 
tity 7 = |ne( 3 )|, where the components of n are portions 
of plaquettes of various orientations and the components 
of e^ 3 ) are cubic roots of unity. That is n = (nx, fly, n>z), 
n x + n Y + n z = 1 and e^ 3 ) = (1, e i27r / 3 , e l47r / 3 ). In the 
case of orientational order, (7) would approach a finite 
value as size of the system increases. In our case it scales 
to instead. 

More insight can be gained from the resonating cube 
correlation functions. In Fig. [5] we plot the function 
G2(r, 0). Unlike the function Gi(r, 0), it clearly captures 
strong persistent correlations in the direction perpendicu- 
lar to the plaquette as well. Our findings for G\, G2 and 
7 can be combined in a coherent picture in which the 
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FIG. 6: Pair correlation function G2(r, 0) for L = 64. 




FIG. 7: (color online) Ordering pattern. Plaquettes form 
resonating cubes which are ordered on a cubic lattice with 
doubled lattice spacing. Different colors represent plaquettes 
of three different orientations. The plaquette orientations are 
uncorrelated over large distances. 



resonating cubes form a lattice, breaking translational 
symmetry, but the orientations of plaquettes on different 
cubes are not correlated over long distances - thus ori- 
cntational symmetry is not broken. This is illustrated in 

Fig.m 

The presence of long range order (LRO) suggested by 
Fig.[S]is supported by a scaling analysis. In Fig.[5]we plot 
differences AG 2 (0,L/2) = G 2 (0,L/2 + I) - G 2 (0,L/2) 
and AG 2 (£/2,0) = G 2 (L/2 + l,0) - G 2 (L/2,0) for L = 
2 P . We see that AG 2 approaches a constant already at 
p ~ 5 — 6. [For larger p, we lose ergodicity, as we show in 
Appendix[D] where we evaluate the relevant equilibration 
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FIG. 8: The curves show scaling of G2 with the size. 
AG 2 (0,L/2) (circles) and AG 2 (i/2,0) (squares) are plotted 
vs log 2 (L). Lines are guides to the eye. 



time r using a binning procedure.] 



C. Origin of the ordered pattern 

Our MC simulations revealed a crystalline phase in the 
classical system of plaquette coverings of 3D square lat- 
tice. The system orders in a cubic lattice pattern with the 
periodicity equal to two lattice constants of the original 
lattice. 

This ordering is tenuous, as evidenced by the small 
value of the order parameter. In fact, when one judges a 
snapshot of the Monte Carlo simulations by eye (Fig. fTO)) . 
no ordering is evident. 

One way to understand the appearance of the long 
range order is to consider a more general model in which 
nrf-dimensional hypercubes are arranged on a n s dimen- 
sional lattice. For RSVP, nd — n s — 3. 

Let us start with n s = 1, hypercubes arranged on a 
chain, Fig. [5] One can then explicitly compute (see de- 
tails in the Appendix [5} the correlation length, which 
grows as y/nd with nd- This result can be visualised by 
noting that a resonating cube can be broken up into two 
domain walls by fixing the plaquettes to be perpendicular 
to the chain direction. This costs an entropy Sd = lnn^. 
The domain walls can then be separated at will without 
disturbing other cubes. In a system of length L, this leads 
to an entropy of S w = ln(„). The correlation length is 
then determined by the smallest L at which the intro- 
duction of domain walls becomes favourable (S w > Sd), 
which happens at L = £ ~ \frid. 

While there is no true long range order in one dimen- 
sion, this result shows that the ordering tendency in- 
creases with rid- In this spirit, the case rid — n s — 2, 
square dimers, being critical, is indeed already closer to 
ordering. Zeng and Henley, for the honeycomb lattice, 
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FIG. 9: (color online) rid dimensional 'hypercubes' resonating 
on a n s — 1 dimensional chain. 



have developed a prescription very similar to the intro- 
duction of a limit of large rid for n s — 2, and they found 
numerically the existence of an ordering transition for 
sufficiently large rid^ 

Thence, whereas for n s = 1, rid — ► 00 yields a critical 
point, n s = 2 requires only a finite to produce order- 
ing. Our results adds the insight that for n s = 3, rid = n s 
is already sufficient to produce ordering. 

Some further understanding can be gained by consid- 
ering higher spatial dimensionality as chains of d = rid- 
dimcnsional hypercubes, transversally coupled in d — 1 
transverse directions. If a chain and its neighbour (in 
the ii direction, say) are in registry, a pair of neighbour- 
ing resonating cubes can gain further free energy when 
their plaquettes are perpendicular to e, . This happens in 
1 out of d configurations for each cube, and the number of 
additional states of the two neighbouring plaquettes be- 
longing to different chains resonating together is d. This 
thus yields an attractive interaction of strength propor- 
tional to 1/d. 

For a d dimensional interaction of strength J x 1/d, the 
critical coupling J would be 0(1). Our interaction, how- 
ever, is amplified by the length of the domains, £ ~ \fd, 
over which individual chains are ordered even in the ab- 
sence of transverse coupling. This means that for large 
d, the transversally coupled chains will order. This es- 
tablishes at least local stability of the plaquette ordering 
for large d, to which our case d — 3 (but not yet the case 
d = 2) is connected. 

Note that this ordering pattern does not imply that, 
even at large d, the plaquettes spend almost all their time 
resonating with their partner on one given cube. (Such 
cubes define the 'site' of the ordering lattice depicted in 
Fig. [?])■ Rather, they spend half of their time resonating 
on the d bonds emanating from sites of that lattice, as 
can be seen from an appropriate adaptation of the above 
argument. 

Let po (pi) be the respective probabilities of finding 
a crystal site (bond) to be occupied by two parallel pla- 
quettes. The probability of a resonating cube to appear 
on any given bond is proportional to p\/d 2 . It should 
be equal to the probability for a resonating cube to leave 
a bond, which is proportional to pi/d. In the large d 
limit other processes are of higher order in 1 j d. Indeed, 
one can evaluate probability p t for a resonating cube to 
appear in locations which are i manhattan steps away 
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FIG. 10: (color online) Snapshot of a plaquette configuration 
on a cubic lattice 

from crystal sites (manhattan distance here is measured 
in original lattice spacings, so a bond is one manhattan 
step away from crystal sites). The relation between pi 
and Pi+i is similar to the relation between po and p±, 
that is p 2 /d 2 ~ Pi+i/d (as long as i <C d). Therefore 
Pi ~ 1/d 2 , while for the number ni of locations at 
manhattan distance i we have rii oc d % . It follows that 
only terms with i = 0, 1 will contribute in d — > 00 limit, 
so a simple relationship holds: p\d — (1 —po)/2. From 
there we find: pq — 1/2 and p\ = l/(4d). 

Because a plaquette ia as likely to be found resonat- 
ing on a bond as it is on a cube, a random snapshot 
of the cubic lattice will not appear to be ordered, simi- 
larly to a snapshot of the face centered cubic (FCC) lat- 
tice (Figs llOlTTI and Sect. IIIip . A proper analysis reveals, 
however, that plaquette correlations on the FCC lattice 
are short range, so the presence of ordering is not a fea- 
ture of resonating plaquettes in general, but should be 
attributed to the lack of frustration on the cubic lattice. 



D. Phase diagram of the RK model 

The above has established that the equal amplitude 
sum, which has zero energy and thus saturates a lower 
bound on the ground state energy at the RK point, is 
crystalline. The QPM Hamiltonian is not ergodic in the 
full Hilbert space — it moves in topological sectors which 
we will discuss elsewhere^ — and thus will properly ex- 
hibit multiple ground states. We can conclude from the 
simulation of the classical ensemble though, that none of 
these ground states will be liquid. 



7 




FIG. 11: (color online) Snapshot of a palquette configuration 
on a FCC lattice. Eight colors are used to represent the dif- 
ferent plaquettes, as defined by faces of an octahedron (see 

FigE). 

The remainder of the phase diagram of the RK model is 
completed in the standard way. For v/t > 2, any configu- 
ration without resonating cubes becomes a ground state. 
The ordered plaquette phase takes over in a first-order 
transition at the RK point v/t = 2. This is unlike the 
dimer model, in which there is a moderately exotic decon- 
fined multicritical transition on the square lattice into a 
solid and into a Coulomb phase on the cubic lattice£&i& 

For — v/t 1, one finds the columnar solid phase, 
which is unique (up to global symmetry operations) in 
this case. Whether there are intervening phases between 
columnar and plaquette solid has not been studied. 



III. N=3 QPM ON THE FACE-CENTRED 
CUBIC LATTICE 

To see if the absence of disordered phases is a general 
feature of such plaquette models, we have also investi- 
gated the case of the face-centred cubic lattice. There, 
the elementary plaquettes are equilateral triangles, the 
corners of which are defined by a triplet of mutually 
nearest- neighbour sites. The QPM now includes reso- 
nance of a pair of plaquettes the six vertices of which 
define an octahedron shown in Fig |121 

In Fig. 021 we plot the plaquette correlations of the 
classical model corresponding to the RK point for var- 
ious inequivalent combinations of pair orientations and 
separations in the [Oil] direction. There are six inequiv- 
alent correlation functions for four types of plaquettes 




FIG. 12: (color online) Octahedron defining the plaquette 
types on the FCC lattice. 

adjacent at an octahedron vertex. All correlators decay 
with distance to the value given by uncorrelated plaque- 
ttes. This decay is a rapid exponential, with correlation 
lengths of less than a lattice spacing, ranging from ~ 0.3 
(diagonal component) to ~ 0.7 (off-diagonal component) 
of nearest neighbor distances. 

This demonstrates that a classical hard-core plaquette 
model is in a strongly disordered phase and we would 
be tempted to conclude, by standard arguments, that we 
have identified a point in an RSVP phase that sits at a 
first order transition to a staggered crystal. The connec- 
tion of the classical sum to the RK point of the QPM, 
however, is more complicated than in the case of the cubic 
lattice. The plaquette pair flip is very constrained and by 
itself breaks up the Hilbert space into a very large number 
of largely frozen states with local quantum fluctuations 
so that it is only the sum of these individual states that 
exhibits liquidity. It seems very likely that the inclusion 
of resonance dynamics on larger clusters will remove this 
pathology and establish the RSVP state. However, at 
present we are unable to explicitly solve this problem. 

IV. CONCLUDING REMARKS 

Can the above quantum plaquette models be realised 
in experiment? Conceivably, they can arise in the case of 
low energy Ising degrees of freedom which reside on the 
midpoints of plaquettes of a lattice, and which obey the 
appropriate hardcore constraint. Their simplest quan- 
tum dynamics would then be that of a QPM. 

Regarding our original motivation of finding degrees of 
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the same lattice i 39 i Similarly, attractive interactions in 
the simple cubic dimer model lead to an ordered phase, 
albeit one with columnar order~ rather than the pla- 
quette order observed here. 

Finally, our work on the RK point also provides an 
entry to the study of classical hard plaquettes in d > 2. 
In particular, the topological and ergodicity properties 
of these models should present an interesting topic for 
further research. 
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freedom with a native SU(iV) symmetry, we first observe 
that following initial work on dimer ground state ensem- 
bles of Klein models 3 ^, there now exist prescriptions for 
local SU(2) invariant Hamiltonians which produce RK 
quantum dimer models in a controlled fashioni 17 i 18 i 38 

However, no simple generalisation of Klein models for 
SU(4) models with plaquette ground states is known yet, 
in part for the reason that was pointed out in Ref. [36| 
for dimers on the square lattice. There, the projectors 
occurring in the Hamiltonian allow ground states con- 
taining not only dimers on nearest neighbour but also 
on next nearest neighbour bonds. In an anlogous man- 
ner, for SU(4) plaquette ground states, we are also lead 
to other (e.g. non-planar) units of four spins forming a 
singlet. Actually realising an SU(4) degree of freedom 
is of course also no trivial matter in d > 1, as combin- 
ing an orbital pseudospin 1/2 degree of freedom with an 
SU(2) to form an SU(4) invariant Hamiltonian requires a 
fine-tuning of interaction parameters which is not imme- 
diately suggested by the natural symmetries of orbitals 
making up the pseudospin 1/2. Similarly, whether or not 
there exist lattice models realising an SU(3) QPM on the 
face-centred cubic lattice that exhibits an RSVP phase 
is left for further investigation. 

It is worth noting that SU(4) spins do of course not 
present the exclusive entry point to plaquette degrees of 
freedom. We have already referred to work on spin-3/2 
SU(2) systems that give rise to plaquette o 29 i 30 . One can 
also consider that in the same sense in which dimers can 
be thought of as low-energy composite degrees of free- 
dom consisting of a pair of resonating SU(2) spins with 
S = 1/2, a plaquette might be considered as a compos- 
ite of a pair of resonating dimers. There is thus nothing 
fundamentally untoward about resonating valence pla- 
quettes - this view being reinforced by the findings of 
Ref. 1261 who found plaquette ordering with the same 
VT~2 x \/l2 supercell proposed for resonating dimers on 



APPENDIX A: PLAQUETTES ON THE CHAIN 

We use the transfer matrix method to compute cor- 
relation function of d dimensional plaquettes on a chain. 
Each plaquette can assume d orientational states: one or- 
thogonal to the direction of the chain and d — 1 parallel. 
We write the Hamiltonian as H — J2i HiA+i an d com- 
pute the transfer matrix T" Q = (a\ exp (— (3Hi,i + i)\a'), 
where \a) represents local degrees of freedom. Provided 
the Hamiltonian is translationally invariant the correla- 
tion function C(x) = (MiMi +x ) of some operator M can 
be written as: 



C(x) 



TtMT x MT 
TrT^ 



N-x 



(Al) 



where N is the length of the chain (periodic boundary 
conditions are assumed). In the thermodynamic limit 
N — > oo we have: 



C(x) =J2(MM\\ a )(K\M\\ ) 



A 



(A2) 



where \X a ) are eigenstates of the transfer matrix and Ao is 
its largest (in absolute value) eigenvalue. The asymptotic 
form of the connected part of the correlation function is 
determined by the next largest eigenvalue A such that 
(A |M|A)(A|M|A ) ^ 0, so that the correlation length 
is given by £ = 1/ In |Ao/A|. In our case all states are 
degenerate, so the transfer matrix is simply proportional 
to 2d — 1 dimensional matrix of the following form: 



T 




(A3) 
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where E is (d — 1) x (d — 1) identity matrix, and all the 
matrix elements of A, B and C are equal 1. The matrix 
C is (d — 1) x (d — 1) square matrix, while ^4 and B are 
1 X (d— 1) and (d — 1) x 1 matrices respectively. It is not 
hard to find that T has only two nonzero eigenvalues: 



A = 1/2 ± Vd- 3/4 



(A4) 



from where we find A/A = -1 + 1/Vd + 0(l/d). The 
first term encodes the fact that the connected part of 
the correlation function oscillates, while the second term 
determines the correlation length £ = \fd. 



APPENDIX B: OVERLAP MATRIX ELEMENT 

In the fundamental representation four basis states 
(flavors) are used to represent an SU(4) spin. The invari- 
ance of the SU(4) singlet plaquette wave function under 
rotations in spin space dictates the singlet to be a fully 
antisymmetrized combination of all four flavors on four 
sites: 



Its 



-L £ (-1)*|1234> 



(Bl) 



P{1234} 



where the four flavors {1234} are permuted in all possi- 
ble ways and p is the parity of permutation P. In each 
term there is one flavor per site, that is flavors are per- 
muted among different sites. The prescription for fixing 
overall sign of the wave function will be given below. We 
consider lattice states which are products of singlet wave 
functions corresponding to each plaquette. We need to 
compute the overlap matrix element between two such 
states. It can be visualized with the help of the transition 
graph - the graph which is made of all partially overlap- 
ping plaquettes. Only those have to be included, as the 
fully overlapping plaquettes multiply the matrix clement 
by 1. The total overlap matrix element is a product of 
overlap matrix elements corresponding to disconnected 
parts of the transition graph, therefore it is sufficient to 
consider only a connected transition graph. 

Two plaquettes can be connected either via a com- 
mon edge or a single common site. We call these strong 
and weak connections, respectively. The common edge 
in the strong connection is called connecting edge, and 
the common sites in weak connections are called connect- 
ing sites. If it is possible to connect two plaquettes by a 
path belonging to the transition graph and going through 
strong connections only, then two plaquettes are said to 
belong to the same strongly connected part. Using the 
fact that each site in the transition graph is shared by 
exactly two plaquettes, it is easy to check that either the 
strongly connected part has exactly four weak connec- 
tions attached, or it has no weak connections and we call 
it a simple graph. 

Each plaquette in the transition graph corresponds to 
a certain term 1 1234) of the singlet wave function. The 
overlap is non zero only when sites shared by different 




FIG. 14: Assignment of flavors which respects the constraint 
of no two sites in a plaquette carrying same flavor. 



plaquettes carry identical flavors in each of the two pla- 
quettes. Thus a nonzero overlap is parameterized by as- 
signing a flavor per site in the transition graph, under 
the constraint that no two sites in a plaquette carry the 
same flavor. This is always possible to do by assigning 
flavors to all sites of the lattice in the following way : 
assign 1 to the point (0, 0, 0), 2 to the point (1, 0, 0), 3 to 
the point (0, 1, 0), 4 to the point (0, 0, 1) and then trans- 
late this configuration by (±1,±1,±1) to cover all sites 
of the lattice (see Fig. [T^|) . It is easy to check that any 
plaquette on the lattice contains four different flavors. 

It has to be decided how to choose signs for plaquette 
singlet wave functions to ensure that the overlap matrix 
elements of the initial and final state of a resonance pro- 
cess are positive. It is straightforward to check that it is 
enough to require the term 1 1234) in each singlet wave 
function, whose flavors coincide with our prescription, to 
be positive. 

The assignment of flavors to the sites in the transition 
graph is not unique. This is most easily seen by making 
a connection between this assignment and a loop model, 
which in turn establishes the link to the four-colouring 
model mentioned in the text. We partition the sites in 
the graph (for a given allowed assignment of flavors) in 
the following way. Connect sites with flavors 1 and 2 
which belong to the same plaquette. Do the same with 
the sites carrying flavors 3 and 4. In this way, all sites 
will be divided into noncrossing loops. Notice that in 
each loop two flavors can be exchanged without violating 
the constraint. The configurations of flavors which can be 
connected by a sequence of such transformations are said 
to belong to the same loop decomposition^ The num- 
ber of such configurations is 2 n ' , where n; is the number 
of loops in the decomposition. Some loops contain weak 
connections, while others do not - we call them nonlocal 
and local loops, respectively. From the flavor constraint 
(any plaquette includes all four flavors) it follows that 
a nonlocal loop goes through at least two strongly con- 
nected parts and has at least length / = 4, while a local 
loop belongs to a single strongly connected part and has 



10 



length 1 = 2. Obviously, ni = tin +n n i, where nu and n„i 
is the number of local and nonlocal loops correspondingly. 
Taking into consideration that n n i < n s and n l u < /, — 1, 
where n l u and li are the number of local loops and the 
number of plaquettes in the ith strongly connected part 
(the case n s = 1 is not included here), we can write the 
following inequality for the contribution S to the overlap 
matrix element from a given loop decomposition: 



where rii is the number of vertices which undergo the 
joining procedure having i colored connections, and thus 
is positive. From Eq. IB4I we have 2no + n\ = n^ + 2n^, 
so the maximum in Eq. IB3I is achieved when no and 713 
are minimized, while n\ and 124 are maximized. In the 
best case n\/2 = n^ = n s /3, and M(n s ) = 3 2rt =/ 3 . Since 
n s < I, the bound for S equals the bound for S times 

M(iy. 



2 ni < 6- 



(B2) 



The prefactor of 1/6 accounts for various ways of dis- 
tributing four flavors among two classes of loops. Notice 
that for a simple graph nu = I and n n i = < 1 = n s , 
so the above result still holds in that case, with the in- 
equality turning into an equality. 

In the next step we derive a bound on the number 
of loop decompositions. For the sake of compactness we 
use terms "vertex" and "connection" in place of "strongly 
connected part" and "weak connection" below. Four con- 
nections belonging to the same vertex are joined by non- 
local loops. This can be done in three different ways. 
Therefore there are at most 3 n " loop decompositions. 
Clearly, not all decompositions are allowed, because for 
two nonlocal loops going through the same vertex (there 
are at most two such loops) a choice of flavors of one 
of the loops uniquely determines flavors of the second 
loop. This restriction can be formulated as a problem 
of coloring loops with just two colors, let us say black 
and white. We start with no nonlocal loops defined, that 
is all connections disjoined and not colored. We will be 
joining four connections (belonging to the same vertex) 
at a time, and color them according to the way they are 
joined. Without loss of generality, we choose to do the 
joining procedure on those vertices which already have 
colored connections. 

In this way colors are always defined uniquely (apart 
from the first step). It is easy to check that there are at 
most three ways to join connections when none or one of 
them has been colored. There are at most two ways of 
joining connections when two of them have been colored. 
There is at most one way of joining, when three or all 
connections have been colored. Thus the upper bound 
on the number M of loop decomposition can be written 



{/O\"0+™1 /1 \ "3+114 
2 " s (2) u 

with the following constraints imposed 



(B3) 



^ rij = n s 

i=0 
i=4 

irii = 2n s 



(B4) 



S< 6 



3* 



(B5) 



APPENDIX C: MATRIX ELEMENT 
COMPUTATIONS FOR SU(4) QPM 

We use I a) notations for the orthonormal states and \ a) 
for the original states. We only need to consider states 
which differ by orientation of a single pair of plaquettes, 
we thus name the states according to the orientation of 
the pair in the YZ, ZX and XY planes as \X), \Y) and 
\Z). We first compute matrix elements (a\H\/3). It is 
easy to check that: 



Hijlaiijkl)) = -\a(ijkl)) 



(CI) 



where Hy = E m „^(i)5™(j) is the part of the 
Hamiltonian acting on the bond (ij). The location of 
the plaquette is explicitly shown in \a(ijkl)) and 
are SU(4) generators. 22 In the fcrmionic representation 

S<m — a ln a n: s0 'S'mlA*) = <Wl m ) where P" is the fla- 
vor. The term S"™ (i)S™(j) simply exchanges flavors 
m and n on sites i and j from where the Eq. ljCip fol- 
lows. We note that only Hij with i and j belonging 
to any of the two plaquettes need to be considered. 
When a ^ f3 there is always a plaquette containing both 
i and j, therefore for any bond (a\Hij\/3) = —x and 
(a\H |/3) = —12a; because twelve bonds belong to the pair 
of plaquettes. When a = j3 only eight bonds belong to an 
individual plaquette, these bonds give (a\Hij\a) = — 1. 
For a bond shared by the plaquettes one can estimate 



a \{a\P}\HP\H K 



= -2x l 



So we 



i=0 



{a\H, 

have \a\H\a) = -8 - 8x 2 . 

Let us compute, for example, (X\H\X) and (X\H\Y). 
To the lowest order \X) = \X) - §(|Y) + \Z)) and 
\Y) = \Y) - f {\X) + \Z)). Then (X\H\X) = (X\H\X) - 
%({X\H\Y) + {X\H\Z) +h.c.) + ... = -8 + l6x 2 + 0(x 3 ). 
The constant term —8 is actually independent of pla- 
quette configuration and can thus be dropped. For the 
off diagonal element wc find (X\H\Y) = (X\H\Y) - 
§{(X\H\X) + (Y\H\Y)) + ... = -5x + 0{x 2 ). Thus we ob- 
tain for the kinetic and potential term in the RK Hamil- 
tonian t — 5x and v = 16x 2 . 
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FIG. 15: The standard deviation <t(t) of (G(0, L/2)) T plotted 
vs Log 2 r for L — 48. 



TABLE I: Equilibration time 



L 


8 


16 


32 


64 


T 


2 4 = 16 


2 7 = 128 


2 11 = 2048 


> 2 15 = 32768 



APPENDIX D: DETERMINATION OF THE 
EQUILIBRATION TIME 



however, the measurements are not independent, the er- 
ror cannot be reduced beyond certain limit, as more mea- 
surements do not add any new information. This allows 
to estimate correlation time corresponding to some phys- 
ical quantity O by measuring fluctuation of the average 
value (0} T computed over the time span r. More pre- 
cisely, we measure (0} T over m successive intervals of 
time of length r. From m average values we compute the 
standard deviation <t(t) of the distribution of (0) T . In 
the same run we can measure (0) T /2 in 2m time intervals 
by subdividing each r interval in two. The correspond- 
ing standard deviation a(r/2) is thus computed from 2m 
measurements. The procedure of subdivision (of creating 
smaller bins, hence the term binning) can be continued 
as far as needed. By plotting cr(r) vs Log 2 r we estimate 
the equilibration time r eq as the time where o(r) stops 
changing as r is lowered. Alternatively, one could keep 
the number m of intervals fixed, while decreasing the 
number of measurements in each interval, by doubling 
the timing between successive measurements r. Again, 
the standard deviation cr(r) would depend on r only for 
t > r eq . This is shown in Fig. [T51 where we plotted cr(r) 
for L = 48, for the correlation function G(0, L/2). The 
equilibration time can be estimated to be r eq ~ 2 14 . The 
equilibration times estimated in this way for other sizes 
are summarized in the Table HI 



The error of statistically independent measurement de- 
creases with the number of measurements n as . If, 
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